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, Abstract. The 21 cm emission of neutral hydrogen is the most promising probe of the epoch of reionization 

H ' (EoR). In the next few years, the SKA pathfinders will provide statistical measurements of this signal, such 

as its power spectrum. Within one decade, SKA should produce a full tomography of the signal. Numerical 
simulations predicting these observations are necessary to optimize the design of the instruments and help in 
the interpretation of the data. Simulations are reaching a reasonable level in terms of scale range, but still often 
rely on simplifications to compute the strength of the signal. The main difficulty is the computation of the spin 
temperature of neutral hydrogen which depends on the gas kinetic temperature and on the level of the local 
Lyman-Q flux (The Wouthuysen-Field effect). A Ts > Tcmb assumption is usual. However, this assumption 
does not apply early in the reionization history, or even later in the history as long as the sources of X-rays are 
too weak to heat the intergalactic medium significantly. This work presents the first EoR numerical simulations 
including, beside dynamics and ionizing continuum radiative transfer, a self-consistent treatment of the Ly-a 
. radiative transfer. This allows us to compute the spin temperature more accurately. We use two different box 

' sizes, 20 h~ x Mpc and 100 ft -1 Mpc, and a star source model. Using the redshift dependence of average quantities, 

, maps, and power spectra, we quantify the effect of using different assumptions to compute the spin temperature 

and the influence of the box size. The first effect comes from allowing for a signal in absorption (i.e. not making the 
Ts 2> Tcmb approximation). The magnitude of this effect depends on the amount of heating by hydrodynamic 
shocks and X-rays in the intergalactic medium(IGM). With our source model we have little heating so regions 
seen in absorption survive almost until the end of reionization. The second effects comes from using the real, 
local, Lyman-a flux. This effect is important for an average ionization fraction of less than ~ 10% : it changes 
the overall amplitude of the 21 cm signal, and adds its own fluctuations to the power spectrum. 
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° ' b the Ihomson scattering of LMB photons by free electrons, 
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tion fraction. The second constraint on reionization arises 
from the spectra of high-z quasars. The Gunn-Peterson 
trough appears as a sharp feature in the quasar spectra 
when the neutral fraction in the IGM surrounding the 
quasar rises above 1%. From observations this happens at 
z greater than ~ 6 (Fan et al. 2006 ). The reduced error bar 
on the value of r in the WMAP5 results, together with the 
quasar constraint, strongly favor an extended reionization 
history. 

Much information is coded in the 21 cm emission. First, 
the nature of the sources (Pop II, Pop III stars or X-ray 
sources: QSO, SNe or X-ray binaries) and their formation 
history determine the evolution of the sky-averaged signal 
with redshift. It may also be possible to pinpoint specific 
events in the reionization history like an average ioniza- 
tion fraction of 0.5 using the redshift dependence of the 
power spectrum (Mellema et al. I2006bl Lidz et al. l2007b|) . 
This type of information will be available in the next few 
years from pathfinders like LOFAR or MWA. Within one 
decade, the Square Kilometer Array (SKA) will deliver a 
tomography (maps at many different redshifts) of the sig- 
nal at 1 comoving Mpc resolution. If the different physics 
can be separated in the signal (Barkana & Loeb 2005a, 
McQuinn et al. |2006j) , we will be able to extract infor- 
mation about the growth rate of large scale structures 
directly. If not, it is possible to use them as independent 
source planes for gravitational lensing by lower z clusters 
and still deduce information on the structure growth rate 
(Metcalf & White UDDTJl. 

Numerical simulations of the 21 cm signal are useful 
to explore the different possible source models and de- 
rive constraints for the design of the future instruments in 
terms of frequency range, bandwidth, angular resolution 
and sensitivity. Comparison with observations will then 
enable us to discriminate between source models and for- 
mation histories. Predicting the 21 cm signal requires the 
knowledge of three local quantities : the baryonic density 
field, the ionization fraction of hydrogen and the Ly-a 
flux (to compute the spin temperature). Moreover, high 
resolution is critical to obtain realistic, non-spherical ion- 
ized regions arising from the high photon consumption in 
dense small scale structures, and large simulation boxes 
(> 100 hr 1 Mpc) are necessary to correctly sample the 
large-scale clustering in the source distribution and the 
abundance of rare sources (Barkana & Loeb 120041 Ilicv 
et al. I2006j) . With the development of 3D radiative trans- 
fer codes (Gnedin & Ostriker [19971 Gnedin & Abel [20011 
Nakamoto et al. [20011 Maselli et al. [20031 Razoumov & 
CardallOSMl Mellema et al. I2006al Rijkhorst et al. 120051 
Susa [20061 Whalen & Norman [20061 Trac & Cen [20071 
Pawlik & Schave [20081 Alt ay et al. I2008[) and increases in 
computational power, it has become possible in the last 
few years to predict the 21cm signal. But some compro- 
mises still have been necessary. The most common ap- 
proximation is to assume Ts 3> Tcmb, where T$ is the 
hydrogen spin temperature and Tqmb is the CMB tem- 
perature. In this regime the 21 cm brightness tempera- 
ture is independent of the spin temperature : no infor- 



mation about the gas temperature and local Ly-a flux is 
needed. Obviously, this approximation fails if either the 
coupling of Ts to Tk, the gas kinetic temperature, by 
Ly-a is weak (early and/or far from the sources), then 
Ts ~ TqmBi or if the gas is not significantly heated in 
the voids by X-rays. In other words, this approximation 
removes any possible absorption regime. However, the sig- 
nal seen in absorption has the potential to be stronger 
than in emission (Furlanetto et al. 2006a). Kuhlen et al. 
p006j) . Shapiro et al. ([20Q6J, Gnedin & Shaver pOOIj) and 
Santos et al. (2007J have probed the absorption regime in 
numerical simulations. In Kuhlen et al. and Shapiro et 
al., this is done taking only the role of collisions into ac- 
count in determining the 21 cm signal, i.e. ignoring the 
role of Ly-a pumping or prior to the existence of any Ly- 
a source. Gnedin & Shaver and Santos et al., in turn, add 
the effect of a Ly-a flux in their analysis, but still use a 
drastic simplification : a homogeneous Ly-a flux changing 
with redshift. Considering a simple r~ 2 profile around the 
sources to compute the Ly-a flux fluctuations, Barkana 
& Loeb (|2005b[) have shown that clustering and Poisson 
noise in the source distribution modify the 21 cm power 
spectrum. Furthermore, several recent works (Semelin et 
al. [20071 Chuzhoy & Zheng [2007J Naoz & Barkana HH]) 
have shown that 3D radiative transfer effects modify the 
expected r~ 2 dependence of the flux in a homogeneous 
medium, and even more so in an inhomogeneous density 
field. The present work includes a full modeling of the Ly- 
a radiative transfer and examines how the resulting 21 cm 
signal compares to the one using a homogeneous Ly-a flux. 

Often working with the Ts S> Tcmb approxima- 
tion, the early simulations used box sizes equal to or 
smaller than 20/i _1 Mpc (Ciardi & Madau [20051 Gnedin 
& Shaver [20041 Furlanetto et al. [20041 Salvaterra et al. 
[20051 Valdes et al. [20061) . Recent works put the effort on 
simulating larger boxes, usually ~ 100 h" 1 Mpc, while try- 
ing to preserve a good enough resolution (Mellema et al. 
I2006bl Zahn et al. [2 007J Lid z et al. I2007al McQuinn et 
al. [20071 Lidz et al. I2007bl Iliev et al. l2008j) . For the 
same resolution, radiative transfer simulations are usu- 
ally more demanding than dark matter dynamical simu- 
lations. Consequently, a common strategy is to run very 
high resolution DM simulations, derive the baryonic den- 
sity field with a constant bias, and compute a clumping 
factor within each larger scale radiative transfer cell. This 
is a method to include the effect of minihalos (Ciardi et 
al. [2006) as subgrid physics. In this method, the least ro- 
bust assumption is the use of a constant baryonic to DM 
density ratio : while correct on a large scale, it fails on 
the scale of clusters where hydrodynamics and heating/ 
cooling processes play an important role. In this work, we 
rely on self-consistent hydrodynamic simulations and do 
not use a clumping factor. 

Section 2 presents the radiative transfer code used, 
Section 3 describes the simulation runs, which are ana- 
lyzed in Section 4. Section 5 summarizes the main results. 
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2. The radiative transfer code : LICORICE 

The LICORICE code integrates three main parts, a 
TrecSPH dynamical part which is not used for this work 
(see section 3.1), and two radiative transfer parts for the 
ionizing continuum and the Lyman-a line. All three parts 
are parallelized using OpenMP. 

2.1. Ionizing continuum radiative transfer 

We use a 3D Monte Carlo ray-tracing scheme. The ra- 
diation field is discretized into photon packets which are 
propagated individually on an adaptive grid and inter- 
act with matter. The numerical methods implemented in 
LICORICE are similar to those presented by Maselli et 
al. (|2003l) for the code CRASH. There are, however, a 
number of differences which are described in the following 
subsections. 

The version of LICORICE used in this work does not 
include H2 formation, He ionization nor diffuse radiation 
from recombination. 

2.1.1. Adaptive grid 

LICORICE uses an adaptive grid based on an oct-tree . 
There is a difference between the usual oct-tree and our 
adaptive grid. The oct-tree keeps splitting a mother cell 
into 8 children cells if it contains more than 1 particle. 
However, the grid for RT stops splitting a mother cell if 
it contains less than a tunable parameter, N max particles. 
Therefore, the RT cells can contain from to N max par- 
ticles while the cells in the oct-tree contain either or 1 
particle. 

For equivalent spatial resolution, such an adaptive grid 
requires lower memory and CPU-time than other grid- 
based RT codes. In this paper, we adopted N max = 30 
for the continuum ionization part, which makes the min- 
imum RT cell size of 1 h~ 1 kpc at the end of the S20 
simulation(z=5.6) and 200 ft, _1 kpc for the S100 simula- 
tion(z=6.6) 

2.1.2. SPH-particle density and grid-cell density 

We compute the gas density at each particle's position 
using the SPH smoothing kernel (Monaghan 1992). We 
use on average 50 neighbor particles, N ne i g hbor, to com- 
pute the SPH smoothing kernel. We call it particle den- 
sity and consider that the particle density is close to the 
physical density of the fluid (Hernquist & Katz ll989"j) . We 
define another density: the average density of an RT cell 
that can contain several particles. We call it RT cell den- 
sity and use it to estimate the optical depth of each RT 
cell. When a photon packet travels through an RT cell 
with optical depth At, the absorption probability is given 
by P(At) = 1 — e _Ar . For each RT cell crossed, we de- 
posit a fraction of its initial photon content oc P(At). 
Then, we distribute the absorbed photons and energy to 
each particle in the RT cell proportionally to its HI mass. 



Time interval between two snapshots 
M-eg^ \ 



j J J Radiative transfer time step 

Adaptive update Regular update 

! At c oi! 

-' ■= g j Recombination and 

cooling time step 

Fig. 1. Schematic representation of the different time 
steps used in LICORICE. 



We do not take into account any form of sub-cell density 
fluctuations. Indeed, particles in one RT cell have similar 
particle density because we limited the maximum num- 
ber of particles in one RT cell at N max = 30 which is 
smaller than Nneighbor, the number of neighbor particles 
of the SPH smoothing kernel. Choosing a smaller N max 
requires using more photons, more memory and updat- 
ing cells more frequently. Updatings cells can be espe- 
cially costly around the sources were cells recieve many 
photons. We use the particle density to update physical 
quantities such as the ionization fraction and temperature. 
Therefore, we compute the photo-ionization rate and tem- 
perature individually for each particle even if they are in 
the same RT cell. Recombination and collisional ionization 
are also computed individually for each particle. 

This scheme is especially relevant if the density field is 
not static i.e. an RT cell does not always contain the same 
set of particles. 

2.1.3. Adaptive time integration 

The integration time step for updating the photoion- 
ization and photohcating rates, Aip>T, is adaptive. We 
use a typical A£ h t much smaller than the time inter- 
val Ai snap between two snapshots of the dynamical simu- 
lation. Furthermore, recombination, collisional ionization 
and cooling are treated with an integration time step 
Ai coo i which is much smaller than AiRT- The photoion- 
ization and photoheating rates are kept constant over all 
the Ai coo i in one A£rt ■ The relation between the different 
time steps is summarized in Fig. [T] 

The ionizing continuum radiative transfer is the repe- 
tition of two steps. The emission and propagation of 7V p h 
photon packets during the time interval Ai rcg (for contin- 
uum transfer an infinite velocity of light is assumed) and 
the update of physical quantities. 

— i) Step 1 : All the sources emit a number iV p h of photon 
packets along directions chosen at random. A photon 
packet leaves a fraction of its content of photon number 
and energy in each RT cell it passes through. Each 
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RT cell records the number of absorbed photons and 
energy during the emission of the 7V p h photon packets 
— ii) Step 2 : Then, the ionization state and the temper- 
ature are updated with the integration time step Ai rcg 
for all particles in all RT cells. This update is applied 
even to the RT cells that absorbed no photon during 
Ai reg . This is necessary in order to consider the effect 
of collisional ionization, recombination and adiabatic 
expansion of the gas (which is interpolated between 
snapshots). In view of the 21 cm emission, the change 
of temperature due to the adiabatic expansion of the 
Universe is especially important in diffuse regions, in- 
dependently of the photo-heating rate. 



However, this regular update is not appropriate where 
the flux of photons is too high (around sources for exam- 
ple). In this case, the number of absorbed photons can 
exceed the total number of atoms available for ionization 
in the RT cell, or, in other words, the fixed update time 
step At rcg is much larger than the characteristic time scale 
of photo-ionization. So, whenever the number of absorbed 
photons accumulated in a cell reaches a pre-set limit, for 
example 30% of the total number of neutral atoms in the 
cell, we update the physical quantities of the cell with the 
appropriate integration time step A£rt(< Ai reg ). 

When N p h photon packets have been emitted, the inte- 
gration time is synchronized as described in step 2 , using 
At Teg for the particles which were not updated since last 
synchronization and with the time elapsed since the last 
update for the others. The relative values of the different 
time steps are specified in section 3.2.1. 



2.2. Ly-a line radiative transfer 

Computing Ly-a transfer is necessary to evaluate accu- 
rately the Wouthuysen-Field effect (see section 3.4) in the 
determination of the 21 cm brightness temperature. 

The part of LICORICE which deals with the Ly-a 
line radiative transfer shares some features with its ion- 
izing continuum counterpart : mainly the Monte Carlo 
approach and the adaptive grid. For the Ly-a transfer, 



Nn 



12 is used. This is a smaller value than for the 



ionizing continuum because Ly-a photons have no feed- 
back on the gas: RT cells do not need to be updated at all. 
There is no direct CPU cost connected to using a smaller 

The methods for resonant line scattering were pre- 
sented in Semelin et al. (2007]). The implementation is 
similar to those of other existing codes (Zheng & Miralda- 
EscudeFMH Cantalupo et al 120051 Dijkstra et al 120051 
Verhamme et al. 2006; Tasitsiomi 2006). Necessary im- 
provements have been implemented in our code for the 
simulations presented in this paper. Let us describe them 
briefly. 



2.2.1. Taking the full effect of expansion into account 

Cosmological expansion is usually introduced in Ly-a 
transfer codes by including a Doppler shift associated with 
the Hubble flow velocities and by using a constant value 
for the expansion factor throughout the photon flight. This 
works well for simulation boxes of moderate size (a few 
10 comoving Mpc). However, during reionization, a pho- 
ton emitted just below Ly-/3 will travel several hundreds 
of comoving Mpc before it redshifts to Ly-a. Moreover, 
box sizes larger than 100 comoving Mpc are necessary to 
avoid underestimating the fluctuation power spectrum of 
the 21 cm line emission. Consequently we took the neces- 
sary steps to include the full effect of cosmological expan- 
sion in our Ly-a simulations : 

— We use the comoving frequency. The usual approach is 
to use the frequency in the reference frame of the center 
of the simulation box and convert back and forth to the 
frequency in the gas local rest frame by a Doppler shift 
including peculiar and Hubble flow velocities. Instead 
we use the frequency in the local comoving frame as the 
reference frequency and include only peculiar velocities 
in Doppler shifts. If a photon enters a cell with local 
comoving frequency at time ij n and exits the cell (or 
undergoes scattering) at time t out , its outgoing local 
comoving frequency v out is simply given by the general 
formula : 



«(tjn) . 

a (tout) 



(1) 



where a(t) is the expansion factor of the universe. As 
explained in Semelin et al. (|2007|) . this avoids second 
order errors in — where da is the variation of the ex- 

a 

pansion factor during the flight of the photon through 
the whole simulation box. However, between the in- 
coming and outgoing comoving frequency in a cell, we 
do make a linear interpolation to compute the comov- 
ing frequency along the photon path. This is necessary 
to use Eq. 13 in Semelin et al. (2007) which allows a 
correct computation of the optical depth even when 
the photon redshifts all the way though the Ly-a line 
within the current cell. The point is that in the current 
version, errors do not add up from cell to cell. 
— The second benefit of using this formulation is that we 
track the variation of the expansion factor a(t) during 
the flight of the photon, and we can use the correct 
value at any given time to compute such quantities as 
the gas physical (non-comoving) density. This avoids 
first order errors in ^ when computing the optical 
depth for example. 

In this framework, we naturally account for retarded 
time. The drawback is that we have to propagate over sev- 
eral box sizes (using periodic boundary conditions) some 
photons that were emitted more than 10 snapshots ear- 
lier. This is necessary for the reasons mentioned above 
but also because the source population changes drastically 
over such a time-scale. 
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2.2.2. Specific acceleration scheme 

A photon which redshifts through the Ly-ev line in a hy- 
drogen medium at the universe average density at z ~ 10 
faces an optical depth of r ~ 10 6 (Gunn & Petersonfl.965, 
Loeb & Rybicki 1999). To compute 10 6 scatterings for the 
~ 10 9 photons used in the simulations is impossible. An 
acceleration scheme has to be devised. The usual method 
is the core-skipping scheme (Avery & House T96H1 Ahn et 
al. [2002). A photon entering the core of the line is sys- 
tematically shifted out when it next scatters off a ther- 
mal atom. This scheme reduces the number of scatterings 
and CPU cost by several orders of magnitude while pre- 
serving the shape of the emerging spectrum. However, in 
this work, we do not need to compute the emerging spec- 
trum. All we need is to know where the scatterings occur. 
Consequently we can use an even more drastic acceleration 
scheme. Although the few scatterings that occur in the 
blue wing of the line are important to determine in which 
location the photon enters the core of the line (Semelin et 
al. I2007[) . their contribution to the total number of scat- 
terings is negligible. Consequently, the following method 
is used : we propagate the photon and compute scatter- 
ings until the local comoving frequency is shifted within 
10 thermal linewidths of the line center, in the blue wing, 
at a location Xo- Then we consider that the ~ 10 6 scat- 
terings that this photon is bound to undergo all occur at 
xo. At 10 thermal linewidths off the center, the mean free 
path of the photons is less than 1 kpc and much shorter 
than 1 pc in the core of the line. We can neglect the spatial 
diffusion while the photon remains within the 10 linewidth 
frequency range. As for the scatterings which occur out- 
side the 10 linewidth range, there are just a few hundred 
of them and we can neglect their contribution to the total 
number of scatterings (but not their importance in deter- 
mining where the core scatterings will occur). 

One factor remains to be specified : exactly how this 
~ 10 6 scattering number depends on the local state of the 
IGM. Looking at the optical depth computed for a red- 
shifting photon, three variable quantities can influence the 
number of scatterings : the local neutral gas density, the 
local kinetic temperature of the gas and the current value 
of the Hubble constant. Escaping the line is achieved by 
a combination of frequency shifts from two sources : scat- 
tering off thermal atoms which produce a random walk 
in frequency, and the systematic shift to lower frequency 
between two scatterings due to cosmological expansion. 
How these two effects combine is not clear at first sight : 
although cosmological shifts are systematically redward, 
they are much smaller than the random thermal shifts 
when the photon is still near the core of the line. We ran 
a series of simulations without any acceleration scheme at 
different values of the local parameters, for a single source 
in a homogeneous medium, and we computed the aver- 
age number of scatterings per photon. It turns out that 
changing the temperature in a reasonable range (1 K to 
1000 K) has no influence on the average number of scat- 
terings (variation less than that the statistical error for 



1000 photons). The approach of Loeb & Rybicki CEHM]), 
developed in a K medium, remains valid in the cosmo- 
logical context. This means, as was checked numerically, 
that the average number of scatterings changes as 1/H(z). 
The number of scatterings per photon is obviously propor- 
tional to the overdensity of neutral hydrogen, however, in 
computing the 21cm emission, we are interested in P a , 
the average number of scatterings per atom per second. 
In P a , the contribution of the overdensity cancels out. 
Consequently the following calibrated formula for iV scat , 
the number of scatterings for each photon, was used : 



A^ a .t. = 8 X 10' 



H(z = 10) 
W) 



(2) 



where H(z) is the Hubble constant at redshift z. Many 
non-linear density structures at subgrid levels produce 
deviations from the above formula (Chuzhoy & Shapiro 



3. Description of the simulations 

3.1. Dynamics 

The dynamical simulations have been run using a modified 
version of the parallel TreeSPH code Gadget2 (Springel 
2005) which explicitly conserves both energy and entropy. 

3.1.1. Initial conditions 

The initial conditions are the reference conditions of the 
Horizon ProjectQ. They have been computed using the 
public version Grafic-2 code (Bertschinger 2001]). The cos- 
mological parameters are those of the concordance ACDM 
flat universe based on WMAP3 data alone (Spergel et al. 
I2007j) : n m = 0.24, Q A = 0.76, Q b = 0.042, with a normal- 
ization of the density fluctuation power spectrum given 
by (78 = 0.76. The corresponding Hubble constant at the 
present epoch is Hq = lOOft-kms -1 Mpc -1 , with h = 0.73. 
The initial redshifts (z = 42.7 for the S20 and z = 29.9 
for the S100 simulation) were chosen using the default 
prescription from Grafic-2. It relies on the Zel'dovich ap- 
proximation to predict the growth of the smallest resolved 
structures until they enter the non-linear regime. Then the 
simulations start. Crocce et al. (|2006|) show that this ap- 
proach leaves some spurious transients compared to the 
more accurate second-order Lagrangian perturbation the- 
ory. The corrections are of the order of a few percent on 
such quantities as the density power spectrum of the dark 
matter halo mass function. We believe that, in EoR simu- 
lations, other uncertainties connected to the source mod- 
eling outweigh this issue. 

Using these cosmological parameters, we have run two 
simulations (S20 and S100) covering different comoving 
volumes, respectively 20 3 h~ 3 Mpc 3 and 100 3 h~ 3 Mpc 3 . In 
both simulations, the number of particles is set to 2 x 256 3 
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simulations 


L 


m D M 


^■gas 


£ 




[h- 1 Mpc] 


[ft- 1 M ] 


[h- 1 M Q ] 


[h^kpc] 


S20 


20 


2.6 x 10 r 


5.5 x 10 6 


2 


S100 


100 


3.2 x 10 9 


6.9 x 10 8 


10 



Table 1. Mass and spatial resolution for the two simula- 
tions S20 and S100. 



(~ 33 x 10 6 particles), where half corresponds to gas par- 
ticles and half to dark matter particles. The mass resolu- 
tions of dark matter and gas/stars particles, as well as the 
softening parameter, are given in tab. [T] For comparison 
with other work, the minimum mass of resolved halos is 
- 5 x 1O 1O M for the S100 and - 4 x 10 s M Q for the 
S20 simulations. However, let us emphasize that we do not 
need to identify halos with a halo finder code: sources are 
formed self-consistently in the dynamical code (see section 
3.1.2). 

3.1.2. Gas physics 

Additional gas physics has been added to the public ver- 
sion of Gadget-2. The cooling and heating of the gas are 
computed following the recipe proposed by Katz et al. 
(1996) where the gas is assumed to be optically thin and 
in ionization equilibrium with the UV background (the 
latter corresponds to updated values of the quasar radi- 
ation spectrum computed by Haardt & Madau 1996). In 
addition to the heating and cooling rates given by Katz 
et al. (]1996|), we also have taken into account the in- 
verse Compton cooling using the formula in Theuns et al. 
(|1998[) . These considerations apply to the dynamical simu- 
lation only, not to the radiative transfer simulations. Non 
equilibrium ionization and heating-cooling processes are 
recomputed independently in the radiative transfer simu- 
lations (see section 2.1). 

Star formation is taken into account by transforming 
gas particles into star particles. We have used the classical 
recipe that mimics a Schmidt law : 

d P* _ Pa /o\ 

dt ~ u ' ( j 

where t* is defined by : 

t* = t * ( — J ■ (4) 

Following Springel & Hernquist (2003), we have set 
to* = 2.1 Gyr for the S20 simulation and a physical density 
threshold of p t h = 5.7 x 10~ 25 h 2 g cm -1 . Moreover, it is 
required that the local gas overdensity exceeds 200 in order 
to form stars. 

Supernova feedback is taken into account by simply 
injecting in the surrounding gas particles an amount of 



10 erg per unit solar mass formed. To avoid instanta- 
neous dissipation by radiative cooling (Katz et al. I1996p . 
80% is ejected in kinetic form, and the rest is ejected in 
thermal form. 

3.1.3. Calibration of the star formation history 

To compare the 21 cm signals between the two simulations 
we need very similar star formation histories. However, 
the star formation rate non linearly depends on the nu- 
merical resolution, especially at high redshift (Springel & 
Hernquist 12003) Rasera & Teyssicr 2006). To avoid an ex- 
tremely CPU time consuming fine tuning of to* m the SI 00 
simulation, we forced the star formation rate to follow the 
one obtained in S20. This was performed by adapting at 
each time step the parameter to*, while relaxing the con- 
straint on the minimum density needed to form stars. 

3.1.4. Output frequency 

Radiative transfer simulations are run as a post-treatment 
of the dynamical simulations, interpolating between the 
recorded snapshots. To obtain a reliable interpolation, 125 
snapshots were recorded from z = 41.66 to z — 5.92. 
Conforming to the Horizon Project chosen rule, the inter- 
val between two snapshots is given in terms of the scale 
factor a by : 



3.2. Reionization 
3.2.1. Initial conditions 

The kinetic temperature of the gas is computed in the 
dynamical simulations assuming equilibrium for the ion- 
ization state with a uniform UV background. However, 
it is recomputed in the radiative transfer simulations us- 
ing the photo heating and cooling rate for a non equilib- 
rium ionization state. The radiative transfer for the S20 
(resp. S100) simulation starts with the same initial tem- 
perature as the dynamical simulations, 34. 7K (resp. 17. 3K) 
computed with RECFAST (Seager et al. UM^i . Using a 
uniform temperature is an approximation: from the epoch 
of decoupling (z ~ 150), different regions evolve adiabat- 
ically with different contraction/expansion factors due to 
the growth of density fluctuations. Moreover, the decou- 
pling from the CMB is actually neither instantaneous nor 
homogeneous. Including these temperature fluctuations is 
not a simple task. We can estimate the rms amplitude 
of the fluctuations to be smaller than 10%. We chose to 
ignore them. 

The initial ionization fraction of gas is set equal to zero. 
We used the recombination rate given by Spitzer (jl978) 
and the collisional ionization in Cen (|1992j) . 

Our requirement for choosing the value of the reion- 
ization time step At rcg is that the relative variation of the 
gas temperature anywhere in the simulation box is less 
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Fig. 2. Top panel: Evolution of the fraction of baryonic 
mass locked into stars. Stars are considered young for the 
first 5-20 Myr. Bottom panel: Number of emitted photons 
per hydrogen atom. 
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Fig. 3. Normalized spectral intensity as a function of 
wavelength. 

than 20%. This ensures that the recombination rate does 
not change too fast within A< rcg . To meet this condition 
we used Ai rog = 1/40 Ai snap (resp.1/10 Ai snap ) for S20 
(resp. S100) and At coo i = 1/100 At rog (both for S20 and 
S100). Note that the S20 run required a smaller time step 
since it has a higher resolution. 

3.2.2. Post-treatment of the dynamical simulations: 
the issue of gas heating 

We use the snapshots of the dynamical runs as inputs for 
the radiative transfer simulation. Consequently, there is no 
feedback on the dynamics from the radiative heating of the 



gas. This has two opposite effects. First it slows down to 
some extent the propagation of the ionization fronts since 
the high gas pressure inside the Stromgren sphere is not 
effective. Second the source formation is not quenched by 
radiative feedback : stronger sources speed up the ioniza- 
tion front propagation. The net result of these two effects 
is unclear. In view of other important uncertainties, like 
the source modeling, the feedback of radiative transfer on 
the dynamics is neglected. 

The second consequence of running the radiative trans- 
fer simulation in post-treatment of the dynamical simula- 
tion is that dynamical effects on the temperature of the 
gas due to shock heating are not included in the radiative 
transfer simulation. 

Shock heating is a sensitive issue since it directly af- 
fects the 21 cm signal, possibly transforming absorption 
into emission by raising the gas temperature above the 
CMB temperature. For example, mini-halos with masses 
between 10 4 to 10 8 M Q form very early during the EoR 
and are dense and warm enough from shock heating dur- 
ing virialization to emit 21 cm radiation. Iliev et al. (|2002p 
and Shapiro et al. (|2006|) show that the emission of mini- 
halos can dominate the 21 cm signal prior to the onset of 
Ly-a coupling. Furlanetto and Oh (|2006cp . however, find 
that the contribution of diffuse regions dominates later 
on. The simulations presented here only marginally re- 
solve the most massive of these mini-halos in the best 
case (S20 simulation), and we chose not to try to include 
them as subgrid physics. But we should keep in mind that 
shock heating both below and above our scale resolution 
may affect the 21 cm signal, as long as it occurs in neutral 
regions (mini-halos, filaments). 

The adiabatic cooling of the gas by dynamical expan- 
sion in low density regions is also an important factor re- 
sponsible for the possible strong absorption features in 
the 21 cm maps : we can take it into account easily during 
the post treatment. We use the fact that, below T = 10 4 
K (always true in neutral, 21cm emitting regions), the 
cooling rate of the primordial gas is negligible, and shock 
heating is limited in the low density diffuse regions of the 
IGM. Consequently, we approximate the thermal evolu- 
tion of the gas in the dynamical simulation as adiabatic. 
We compute the variation of internal energy of a particle 
due to expansion or contraction between two consecutive 
snapshots from the variation of the density : 

A it = U2 — u\ = ui ^ (j^j ~ 1J ; (6) 

where u is the internal energy and p is the gas density. 

However, if the variation in density between two snap- 
shots is large, applying it all at once results in numerical 
instabilities. This happens when particles pass through or 
collapse into a halo. To avoid this, the density is interpo- 
lated gradually between two snapshots. 

The other cooling and heating processes are then in- 
cluded self-consistcntly in the radiative transfer Simula- 
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simulations 


L tot [erg/s] 


j esc 


S20 


4.60 x 10 43 


0.08 


S100 


5.75 x 10 45 


0.05 



Table 2. The total ionizing luminosity of a star particle 
and escape fraction of simulations 



tion. We use the cooling- heating rates in Cen (|19921) , and 
a minimum gas temperature of IK. 

3.2.3. Source modeling 

In this work we consider only stellar sources and do not 
study the influence of X-ray sources. The decline of the 
quasar luminosity function at redshift z > 3 suggests that 
stars were the dominant source of ionizing photons during 
the EoR (Shapiro & Giroux fl987|, or Faucher-Giguere et 
al. 120081 for recent observational evidence). However, even 
in small quantities, X-ray sources have the potential to 
affect the 21 cm signal by heating the IGM. Indeed X- 
ray photons can have a mean free path of several 100 Mpc 
through neutral hydrogen during the EoR: they propa- 
gate through the ionization front out to the neutral re- 
gions. Using a homogeneous model Glover & Brand (2003) 
find that a modest number of X-ray sources can raise the 
temperature of the IGM by several tens of K. Moreover, 
Pritchard & Furlanetto (|2007|) show that the heating by X- 
rays is actually not homogenous: it produces fluctuations 
in the gas temperature which leave their own imprint on 
the 21 cm power spectrum. In this work, however, we want 
to focus on the specific influence of Ly-a pumping. X-ray 
sources will be included in a future work. 

All baryon particles have the same mass : 7.6 x 10 6 M Q 
(resp. 9.5 x 10 8 A/ Q ) for the S20 (resp. S100) simulation. 
Accordingly, one star particle corresponds in fact to a star 
cluster or a dwarf galaxy so an IMF has to be assumed. 
We choose a Salpeter IMF, with masses in the range 1M© 
- IOOMq. Then we compute the total luminosity and spec- 
tral energy distribution (SED) for one star particle by in- 
tegrating over the IMF. We use the table in Aller (|1982|> 
to obtain the radius and effective temperature of stars as 
a function of mass. The total ionizing luminosity (tab. [5]) 
of a star particle yields one or two photons per hydro- 
gen atom at z=6.5 (Fig. [2]). The resulting SED is plotted 
in Fig. [3l It can be seen that it differs very little from 
a blackbody spectrum with 5 x 10 4 K effective tempera- 
ture. Therefore our source model is roughly intermediate 
between PopII and PopIII stars. 

3.2.4. Calibrating the global ionization history 

Before making a global calibration using the escape frac- 
tion, we need to consider the effect of the time evolution of 
the stellar population within one source particle. Massive 



stars contribute the most to the ionizing luminosity of 
the source, but live for a shorter time. We use a simple 
model : when a new source particle is created, we assign it 
a random lifetime between 5Myr and 20Myr as a source of 
ionizing photons. The randomization is useful because we 
obtain new sources only with each new dynamical snap- 
shot; a fixed lifetime for all sources and a discreetly evolv- 
ing number of sources with each new snapshot would pro- 
duce distinct steps in the average ionizing flux. Globally, 
the average lifetime of the source is degenerate with the 
escape fraction as far as ionization history calibration is 
concerned, however, it allows us to account for the signa- 
ture of local starbursts (in time and space) followed by 
quiescence. We plot the effect of considering a finite life- 
time for the ionizing source on the star fraction history in 
Fig. [21 It can be seen that the young, ionizing star frac- 
tion fluctuates more strongly at high redshifts than the 
global star fraction. A combination of local starbursts fol- 
lowed by star formation quenching and insufficient mass 
resolution may be responsible for this effect. 

By trial and error, we calibrated the escape fraction 
(tab. [2]) so that reionization is complete, at about z ~ 6. 

3.3. Lyman-a transfer 

We run the Lyman-a simulations as a further post- 
treatment of the reionization simulations. The time inter- 
val between two snapshots is typically 10 Myr at z ~ 9, 
during which the ionization fronts travel a few hundred 
comoving kpc at most. We consider that this time inter- 
val is small enough to permit a de facto interpolation of 
the quantities relevant for Lyman-a transfer : gas density 
and ionization state. 

Potentially, Lyman series photons can heat up the gas 
(see references in Semelin et al. I2007P . But heating occurs 
only at low temperatures (~ 10 K), where the gas pressure 
is insignificant. Consequently, the feedback of Lyman-a 
photons on the dynamics is negligible. However, in models 
with soft source spectra, where the efficient heating of low 
density neutral regions by X-ray photons is absent, one 
consequence of including the heating by Lyman-a photons 
is to increase moderately the temperature in the coolest 
regions. We implemented a local Lyman-a heating using 
the following formula, derived from Eq. 47 and Fig. 3 in 
Furlanetto & Pritchard (|2006b[) : 

dT K 0.8 x a 10 / 0.7 \ 

-dT = ^s- a TT-z Hiz) Wlth Sq ~1 1_ ^J (7) 

We checked that when the 21 cm signal is in absorption 
with an intensity in the 300 mK range, i. e. where the gas 
temperature is a few K, this small amount of heating can 
reduce the signal intensity by up to 100 mK. 

Photons are emitted at frequencies between Ly-a and 
Ly-/3. The spectrum and luminosity of the sources are 
computed with the method described in section 2.2.3 for 
the ionizing continuum. In principle higher Lyman se- 
ries photons also contribute. Indeed as soon as a photon 
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emitted in the continuum hits an atom in a resonant up- 
per Lyman line, it cascades locally down the levels and 
is transformed into a local Ly-a photon. These higher 
Lyman line photons are more concentrated around the 
sources since they need a smaller redshift before they res- 
onate with a line (see e.g. Naoz & Barkana 2008;). Overall, 
including the upper Lyman line photons would increase 
the number of Ly-a photons by a factor of less than 2. 
From the results of the simulations we estimate that this 
could shift the history of the average Ly-a coupling by 
Az — 0.3. It would also add some limited amount of power 
in the fluctuations. In this work we include only the Ly-a 
line. Upper Lyman series lines will be included in a future 
work. 

In the simulations, a fiducial number of = 10 7 pho- 
tons is used between two dynamical snapshots. As the time 
interval between snapshots and the total source intensity 
varies, the physical content of the photons also varies. The 
value of rtph = 10 7 , dictated by the CPU cost of the simu- 
lations, is a lower limit for statistical significance : indeed, 
with our current grid resolution of less than 12 particles 
per cell, about 10% of the cells do not receive any pho- 
tons between two snapshots. These cells are small and are 
mostly located in moderate density filaments and low flux 
environments, far from the sources. They are properly av- 
eraged with some of their neighbors when the fixed-grid 
maps are computed. To assess the impact of this limited 
sampling we also ran a simulation with 6 x 10 7 photons be- 
tween two snapshots down to z = 9.5 for the 20 h _1 Mpc 
box. In this case the fraction of unsampled cells drops to 
~ 1%. It is only when maps are drawn with thin slices, 
~ 200 ft. -1 kpc, that a difference can be seen. 

3.4. Physics of the 21 cm emission 

To compute the differential brightness temperature of the 
21 cm emission in the optically thin limit of the 21 cm line, 
we use the formula : 



ST b = 28.1 mK x m (1 - 
Qb h 



S) 



0.042 0.73 



1 + z 
10 

0.24 



l-Y p 
1 - 0.248 



(8) 



where S is the local overdensity at redshift z, xm the 
neutral hydrogen fraction, Ts the neutral hydrogen spin 
temperature, Tomb the CMB radiation blackbody tem- 
perature at redshift z and f2b, fi TO , h and Y p are the 
usual cosmological parameters. In this formula, the con- 
tribution from the proper velocity gradients have been ne- 
glected (Barkana & Loeb I2005a|) . Mellema et al. (|2006b|) 
include this effect in computing 21 cm line of sight spectra. 
Analytic (McQuinn et al. I2006P and semi- numerical mod- 
els (Mesinger & Furlanetto I2007j) find that the contribu- 
tion of proper velocity gradients to the 21 cm power spec- 
trum is important only at (xhi) < 0.5 on scales smaller 
than the effective ionization bubble size. Since previous 



S20: mass average 
S20: volume average 
SlOO: mass average 
SlOO: volume average 
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Fig. 4. Mass weighted and volume weighted averaged ion- 
ization fraction. 



simulations did not include a proper treatment of the Ly- 
a pumping which is especially important during the early 
phase of reionization, it was not relevant to include the ef- 
fect of proper velocity gradients. In this work, we include 
a full modeling of the Ly-a pumping. In a separate paper, 
the effect of proper velocity gradients will be studied. 

The spin temperature of hydrogen Ts is coupled to 
the CMB temperature Tcmb by Thomson scattering of 
CMB photons, and to the kinetic temperature of the gas 
Tk by collisions and Ly-a pumping. As a result, the spin 
temperature can be written (Furlanetto et al, 2006a) : 



Tr 



T-l _ t CMB j~ x cT K 



x a T, 



c 



with Tn ~ T-. 



K 



(9) 



and 



4P Q T, 



27v4i Tcmb 



(C H + C p + C n ), (10) 



where T* = = 0.068 K {v = 1420MHz), A w = 
2.85 x 10 _15 s _1 is the spontaneous emission factor of the 
21 cm transition, P a is the number of scatterings of Ly-a 
photons per atom per second, and Ch, C p and C e are the 
de-excitation rates due to collisions with neutral atoms, 
protons and electrons. For the de-excitation rates we use 
the fitting formula given by Liszt (pOOT]) and Kuhlen et al. 
(2006). The coupling coefficients x a and x c are computed 
for each gas particle from the simulation data. The final 
output is the value of the brightness temperature for each 
particle in each snapshot. 

4. Analysis 

We will focus our analysis on two aspects. 

First, by comparing our two simulations, we will eval- 
uate the effect of numerical resolution and box size. 
Numerical resolution truncates the physics at small scales : 
this mainly modifies the reionization history and geome- 
try by removing small, dense structures with high recom- 
bination rates. The limited box size, on the other hand, 
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Fig. 5. Evolution of the volume averaged Lyman-a and 
collisional coupling coefficients as functions of the average 
ionization fraction for both simulations (see main text for 
a definition of the coefficients). 



truncates large scales and dampens the 21 cm brightness 
temperature power spectrum on corresponding wavenum- 
bers. 

The second axis of our analysis is to measure the effect 
of including an accurate treatment of the Wouthuysen- 
Field effect in computing the brightness temperature. 

4.1. Evolution of sky-averaged quantities 
4.1.1. Ionization fraction 

The simplest way to characterize the global history of 
reionization is to compute the average ionization fraction 
as a function of redshift. The mass weighted ((£ff)mass) 
and volume weighted ((xh)voi) average ionization frac- 
tions are shown in Fig. 3] for both simulations. For the 
volume (resp. mass) average we use the volume occupied 
by each particle derived from the SPH smoothing length 
(resp. the mass of each particle) as weights in the averag- 
ing procedure. The main feature of this plot is that the 
mass weighted averages are similar for the two simulations 
while the volume weighted averages differ. Indeed, the star 
formation histories of the two simulations were calibrated 
to produce the same number of ionizing photons at any 
redshift. As long as recombination is negligible, this should 
produce identical mass weighted ionization fractions. It 
can be checked in Fig. [J] that below redshift ~ 8 recom- 
bination becomes efficient in the S20 simulation, which 
contains denser small-scale structures due to higher res- 
olution, and the mass weighted ionization fraction drops 
below its S100 counterpart. 

The volume weighted ionization fractions differ at all 
redshifts between the two simulations, because highly ion- 
ized regions of identical mass, close to the sources, are 
denser and occupy a smaller volume in the more resolved 
S20 simulation. McQuinn et al. ([2006) advocate the use 
of the volume-averaged ionization fraction as an evolu- 



Fig. 6. Evolution of the neutral hydrogen spin tempera- 
ture as a function of the average ionization fraction. The 
average is computed either directly from the particle spin 
temperatures or from the volume averaged values of Tk 
and x a . The CMB temperature and volume averaged ki- 
netic temperature are plotted for comparison. 



tion variable instead of redshift. Indeed, depending on ar- 
bitrary choices (all consistent with current observational 
constraints) for the source formation history and source 
nature, reionization can begin at z = 12, 15 or 20 and 
proceed differently. However, for a fixed average ionization 
fraction, McQuinn et al. ( 2006) show that the brightness 
temperature power spectrum depends weakly on the red- 
shift at which this ionization fraction is achieved. This is 
true if a modification of the source formation efficiency is 
the main cause of the change in redshift for a fixed ioniza- 
tion fraction. If the change is caused by using a different 
source spectrum, or modifying the gas sub-grid physics 
(e.g. an inhomogenous clumping factor), this may not be 
true anymore. We believe that the averaged ionization fac- 
tor is indeed a more relevant measure of the evolution 
of reionization than redshift and we adopt their point of 
view. However, instead of the volume averaged ionization 
fraction, which is well suited to comparing different models 
with the same underlying gas density field, we will use the 
mass averaged ionization fraction which is better suited 
for the comparison of simulations with different mass res- 
olution. While the mass averaged ionization fraction is a 
suitable quantity to track the progression of reionization, 
a volume average is more relevant to compute such ob- 
servable quantities as the 21 cm brightness temperature. 
Consequently, all quantities other than the ionization frac- 
tion will be volume averaged. 

We also computed the Thomson optical depth from 
the two simulations ignoring the presence of helium. The 
values are r = 0.060 and t = 0.063 for S20 and S100 sim- 
ulations respectively, i.e. within la of the WMAP3 value 
t = 0.0089 ± 0.030 (see Spergel et al. l2007|) . 
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4.1.2. Spin temperature coupling coefficients 

The spin temperature of neutral hydrogen, Ts, is defined 
by Eq. 7 and depends on two coupling coefficients, x a 
for the Wouthuysen-Field coupling and x% for collisional 
coupling. Until now, x a has never been computed consis- 
tently in simulations of the 21 cm emission. While x a itself 
contains fluctuations which contribute to the power spec- 
trum of the 21cm signal, we study first the evolution of 
the average value of the coefficients which tells us at which 
stage in the reionization history the spin temperature can 
be considered fully coupled to the kinetic temperature of 
the gas. The evolution of the volume averaged coefficients 
{%a)vo\ and (xk)vo\ is shown in Fig. as a function of 

y^H /mass- 

The first conclusion from Fig. [5] is that Xk is 
negligible compared to x a except in the very early phase 
((£ff)mass ~ 10~ 4 ) for the S20 simulation. Logically, the 
S20 simulation, producing denser structures, has a higher 
Xk (indeed Ch, C p and C e in Eq. 8 depend on the density). 
With a better mass resolution and even denser structures 
than in the S20 simulation, Xk would still increase, and 
the ionization fraction when x a starts to dominate would 
also increase. The second piece of information that can be 
derived is the magnitude of the error made by applying 
the usual assumption x a = +00. When x a > 10, the er- 
ror made on 5Tb by assuming x a — +00 is smaller than 
— 10%. This corresponds to (x H ) mass > 0.02 for the S20 
simulation and (x^/) mass > 0.04 for the S100 simulation. It 
may appear that the full coupling occurs very early on. Let 
us emphasize, however, that in this early phase there is a 
strong absorption signal (with our choice of source model) 
compared with the weaker emission signal observed later 
on. In a source model with a harder spectrum (to be stud- 
ied in a future paper) we would have a weaker absorption 
signal. If the sources are still stars ( e.g. using a more top- 
heavy IMF) the full coupling would also occur at a larger 
average ionization fraction since the ratio of ionizing to 
Lyman-alpha photons would be larger. 

4.1.3. Spin temperature 

Fig. [fj] displays the average spin temperature for the S20 
simulation, with the average computed in two different 
ways. First the spin temperature is computed using the 
volume averaged quantities (2fc) vo i, (x a ) vo \ (and (xk) vo \), 
following Furlanetto et al. (2006a) and Gnedin & Shaver 
(|2"0M)) . This quantity_is denoted T s (f K ,x a ). Then the 
more relevant value, Ts, is computed, as the direct vol- 
ume average on the particles of the spin temperature, i.e. 
Ts = (2s) vol. Fig. [6] shows that there is a large differ- 
ence between the two definitions. This is mainly due to 
the fact that the spin temperature is computed through 
the weighted average of inverse temperatures (Eq. 7) . This 
large difference propagates to the average brightness tem- 
perature. Consequently we strongly advocate replacing 
the Furlanetto et al. average traditionally used to define 
the emission and absorption regimes by the more rele- 
vant direct average of the spin/brightness temperature. 
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Fig. 7. The evolution of the volume averaged spin tem- 
perature, kinetic temperature, and CMB temperature is 
plotted as a function of the average ionization fraction. 
The quantities are plotted for both box sizes. 
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Fig. 8. Evolution of the average differential brightness 
temperature for the S100 simulation. The average is com- 
puted in three different ways (see main text and legend). 
The blue curve is the most relevant direct average of the 
particle brightness temperature. 



Of course this is easily done only for simulations. Let us 
mention that with the direct average, the spin tempera- 
ture remains below the CMB temperature at lower z. 

Fig. [7] compares the direct average of the particle spin 
temperature in the two simulations : the results are quite 
similar. The main difference is that the spin temperature 
starts to decouple from the CMB temperature later in the 
SI 00 simulation. Indeed, we use periodic boundary con- 
ditions in the simulations. Thus the Lyman-a sources are 
more homogeneously distributed in the 20/i _1 Mpc box, 
which lacks large scale modes, than in the 100/i _1 Mpc 
box (Barkana & Loeb I2005bj) . In the lOO/i^Mpc box, 
large void regions receive little Lyman-a flux. There, the 
spin temperature remains coupled to the CMB at lower 
rcdshifts. 
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Fig. 9. Comparison of the evolution of the average differ- 
ential brightness temperature for the two box sizes. The 
average is computed in two different ways (see main text). 



4.1.4. Differential brightness temperature 

The differential brightness temperature is plotted in Fig. [8] 
for the SI 00 simulation. The average is computed in three 
different ways : 6Tb(xu, P, Ts) with Ts = Tk for each par- 
ticle, STi,(xh, P,Ts) computing Ts for each particle using 
x a , and STf,, using the same notations as for the spin tem- 
perature. In all these definitions we use a volume weighted 
average. The first definition uses the common Ts = Tk 
assumption : as can be seen in the plot, it fails at early 
times when the coupling is weak. The second definition 
(equivalent to Furlanetto et al. I2006a[) shows an emission 
and an absorption regime. However, with the more real- 
istic third definition, 5Tb, the emission regime disappears 
completely. Indeed, absorption region have a strong sig- 
nal which dominated over the weak, ~ 25 mK signal of 
the emitting regions. Of course, if X-ray sources were in- 
cluded, the cold neutral IGM would be pre-heated and the 
absorption regions contribution would shrink much faster : 
5Tb would show an emission regime and the absorption 
regime would be weaker. Including shock heating would 
go in the same direction although the magnitude of the 
change would be less homogenous and is more difficult to 
estimate. 

To assess the effect of resolution and box size, we plot 
in Fig. [9]the direct 5Tb and 5Tb(xn, p,Ts) for both simu- 
lations. The overall behaviour is similar. In the S20 sim- 
ulation the minimum of 5Tb occurs at a 1% average ion- 
ization fraction compared to a 2% ionization fraction for 
the S100 simulation. This small difference is due in part 
to the slightly smaller coupling of the S100 simulation at 
equal ionization fraction. The values of the minimum 5Tb 
are not very different either. We expect to find resolution 
effects and box size effects in the power spectra. 

We have presented these average values to compare 
with previous studies. However, much information is 
erased when considering average quantities, which makes 



their interpretation of limited scope. We proceed now to 
present maps of the different quantities. 

4.2. Maps 

4.2.1. Ionization fraction 

Maps of the ionization fraction of hydrogen are presented 
in Fig.[TU]for both box sizes at a redshift when (x a ) vo \ = 1, 
which we will call the Wouthuysen-Field coupling redshift 
(hereafter WFCR), and at a redshift when (cc_f/) mass = 0.5, 
which we will call the half reionization redshift (HRR). 
We can see that the WFCR correspond to an early stage 
in the reionization history. The maps are presented for 
two different slice thicknesses : 2 h^ 1 Mpc and 200 h^ 1 kpc. 
The 2 ft, -1 Mpc slice corresponds to a bandwidth 5v ~ 
0.15 MHz at a redshift between 6 and 10, an acceptable 
value for SKA. The 200 h _1 kpc slice is plotted to give a 
feeling of how much power would be added at small scales 
by using narrower bandwidth : indeed, sharper ionization 
fronts are found in this case. Also as expected the S100 
maps show coherent structures on scales up to 50 h" 1 Mpc 
which cannot be accounted for in the S20 maps. 

Usual effect of a lower mass resolution (S100 vs S20) 
is to delay reionization. Indeed small mass halos form first 
and make a large contribution to the ionizing photon bud- 
get at early times (see Iliev et al. I2007|) . In our treatment, 
the delay is removed because we calibrated the source for- 
mation history of the S100 simulation on the S20 simu- 
lation by boosting the formation efficiency. However, this 
procedure can only account for unresolved halos located 
within more massive, resolved, proto-halos. The contribu- 
tion from unresolved halos located in comparatively un- 
derdense regions could modify the morphology of the ion- 
ization maps during the early reionization history. 

4.2.2. Lyman-a coupling 

The influence of including a proper treatment of the 
Wouthuysen-Field effect on the 21cm signal is twofold. 
First, it simply modifies the intensity of the signal, espe- 
cially in the regions seen in absorption. This will be obvi- 
ous in the brightness temperature maps shown in section 

4.2.3. The second effect is to introduce a new source of 
fluctuations in the signal, associated with the fluctuations 
in the local Ly-a flux. Fig. [TT] presents maps of the x a co- 
efficient for both simulations at the WFCR. Around the 
WFCR, the additional fluctuations to the brightness tem- 
perature are expected to be maximal. On a large scale, 
the maps boil down to strong coupling regions around 
sources, superimposed on a uniform weak coupling back- 
ground. The signature of the added brightness tempera- 
ture fluctuations will be based on the typical sizes of the 
strong coupling regions and on the clustering and Poisson 
noise in the source distribution. Keeping this in mind, it 
is interesting to notice that the strong coupling regions 
appear to be much smaller, ~ 2 h^ 1 Mpc, in the S20 simu- 
lation than in the SI 00 simulation where they extend over 
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Fig. 10. Maps of the ionization fraction for the two simulations (right and left) at two different stages of reionization 
(above and below line), the WFCR and the HRR (see text for definition) and for two slice thicknesses (as labeled). 
The color scale is a linear function of the ionization fraction. 
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Fig. 11. Maps of the x a coupling coefficient. We use a logarithmic color scale, and choose the redshift in each simulation 
(S20 and S100) such that (x a ) = 1 (10.25 for S20 and 10.01 for S100). The thickness of the slice is 2 h' 1 Mpc for the 
two panels on the left and 100 h _1 .kpc for the right panel. 



~ 10/i _1 Mpc. To understand this let us remember that 
we apply periodic boundary conditions to the radiative 
transfer code. As a result, any point in the S20 simula- 
tion is surrounded by a rather homogeneous distribution 
of sources, with the closest source at most ~ 10 h^ 1 Mpc 
away. This dampens the fluctuations in the Ly-a flux and 
washes out the outer parts of the strong coupling regions. 
As a result, we expect to find more power at large scales 
in the added brightness temperature fluctuations in the 
S100 simulation than in the S20 simulation. 

The zoom on the source in the S20 simulation shows 
that the apparent spherical symmetry of the strong cou- 
pling regions breaks down at small scales. Semelin et al. 
(2007J) have shown that substantial asymmetry can be ex- 
pected if the density field of the gas is not homogeneous 
(e.g. presence of filaments). Let us notice, however, that 
the slice thickness for the zoom is 100 h _1 kpc, while it is 
2 hr 1 Mpc for the other maps. With such a thin slice the 
noise in the x a evaluation (resulting of the finite number of 
photons used in the Monte Carlo scheme) is not smoothed 
out at all. Any structure outside the main strong coupling 
region is mostly noise (visible on the color version only). 
If structures exist in the x a maps at scales smaller than 
1 hr 1 Mpc and with an amplitude of less than a factor of 3 
above or below the background, we will miss them in the 
noise or smooth them out. From a theoretical point of view 
it would be interesting (but very costly) to produce maps 
with a lower noise level. However, from a practical point 
a view, these scales are below the projected resolution of 
SKA. 

4.2.3. Brightness temperature maps 

The most common (and most drastic) approximation used 
to compute the brightness temperature is Tg Tqmb 
(case 1). In this approximation no absorption region can 
exist. The second approach is to assume full Lyman-a cou- 
pling (x a = oo) (case 2). The third approach, the one used 



in this paper, is to compute self-consistently the local val- 
ues of x a and derive accordingly the local values of T$ and 
Tb (case 3). We plot maps of the differential brightness 
temperature for all three cases and for both simulation 
boxes at the WFCR in the 6 top panels of Fig. [H The 
thickness of the slice is 2 ft, -1 Mpc . At this early redshift, 
when the two first approximations are hardly reliable, the 
differences are striking. Not surprisingly case 1 fails com- 
pletely : it is not tailored to handle absorption regions 
and predicts only emission. More interesting, but also ex- 
pected, case 2 overpredicts the strength of the absorption 
compared to case 3 (300 mK instead of 200 mK) by as- 
suming full Wouthuysen-Field coupling. We used the same 
color scale for all maps to emphasize these differences. The 
strong absorption signal observed in cases 2 and 3 is the 
consequence of our choice of rather soft-spectrum sources. 
The neutral IGM undergoes very little pre-heating by UV 
and X rays. Let us emphasize that the strength of the ab- 
sorption signal is very sensitive to the thermal modeling 
of the gas when Tk <C TqmBj which is not uncommon 
at this redshift with our choice of sources. Indeed, the in- 
cluded Ly-a heating, although very weak (at most a few 
K over the EoR), decreases the strength of the absorption 
signal by 50 to 100 mK. If even a small fraction of hard- 
spectrum sources was included, the absorption strength 
would decrease and turn to emission at later redshifts . 

The six bottom panels of Fig. [T2l shows the same maps 
for the HRR. At this redshift we have (x a ) > 200 for 
both simulations. Consequently the maps for cases 2 and 
3 are almost identical. Comparing with Fig. [10] it can be 
checked that the regions which are still seen in absorp- 
tion are the neutral regions not too close to the ionization 
fronts. Comparing the absorption regions in the S20 and 
in the S100 maps, it can be seen that the signal is some- 
what stronger in the S20 simulation. We have determined 
that a lower gas kinetic temperature in the S20 simula- 
tion, in the voids, is responsible for this effect. Since we 
are talking about temperatures lower than 10 K, a differ- 
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Fig. 12. Differential brightness temperature maps for both box sizes and 3 different hypotheses for computing the 
spin temperature (see main text and description in the figure). The redshift of the 6 top panels is chosen such that 
(x a ) = 1, the redshift of the 6 bottom panels is such that (xu) = 0.5. The temperature color scale is linear, in mK. 
The thickness of the slice is 2 Mpc . 
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Fig. 13. Ratio of the differential brightness temperature 
computed with the fluctuating local value of x a on the 
one hand and the average value of x a on the other hand. 
The color scale is linear. The contour is computed from 
the x a field for the x a = 1 value. 



ence of a few K produce a large effect on the differential 
brightness temperature. The S20 simulation has a higher 
mass resolution which creates higher density contrasts. In 
the voids, the density is lower and the adiabatic expansion 
stronger : the resulting temperature is lower. 

Obviously, in case 1, the absorption regions turn to 
emission. Since emission saturates quickly at a few tens of 
mK, the statistical property of the signal such as the av- 
erage rms fluctuations of the brightness temperature will 
be weaker. At such an advanced stage in the reionization 
history the assumption made in case 1 may actually be 
more realistic: case 2 and 3 would probably look more like 
case 1 if we included X-ray sources. 

Using a homogeneous average value for x a to account 
for the Wouthuysen-field effect ( for example a calibrated 
function of xh) is much easier than computing the full ra- 
diative transfer in the line. So it is important to assess the 
benefit reaped from doing the full treatment. The x a field 
induces its own fluctuation in the brightness temperature 
signal. Fig. [T3l illustrates this point. It shows a map of the 
ratio of the differential brightness temperatures computed 
with the fluctuating local value of x a on the one hand and 
the homogeneous, average value of x a on the other hand. 
This map is computed at the WFCR, for the S100 sim- 
ulation. As expected, the ratio is mostly larger than one 
inside the x a — 1 contour and lower outside. Moreover 
the two differential brightness temperatures differ by up 
to 30%. Consequently it appears that, at these early red- 
shifts when the Lyman-a coupling is not full, it is worth 
doing the full Lyman-a transfer. 



Using the power spectrum as a statistical tool for an- 
alyzing the EoR 21 cm signal is natural because it can 
be directly computed from the visibilities measured by 
the radio interferometers such as LOFAR, MWA or SKA. 
Morales & Hewitt (|2004p have shown that using the full 
3D power spectrum rather than the 2D angular power 
spectrum will improve the sensitivity of the observations 
and make the process of foreground removal easier. 

4.2.4. Strength of the Lyman-a coupling fluctuations 

Fig. [14] shows the power spectrum of the quantity z a = 
j^r- for both simulations at the WFCR. Indeed directly 
plotting the spectrum of x a is not interesting : the x a field 
is dominated by the very strong values near the sources 
and has a power spectrum similar to that of a distribu- 
tion of Dirac functions. However, these very high values of 
x a are not especially interesting : they are just values for 
which the coupling saturates. z a , which is just the weight 
applied to the inverse kinetic temperature in Eq. 9 (ne- 
glecting the contribution of x c ), varies between and 1 
and contains more relevant information about the magni- 
tude of the coupling. As can be seen in Fig. [14] the point 
sources still dominate the small scales. More interesting 
is the increase of the power spectrum toward large scales, 
which is barely visible for the S20 simulation but obvi- 
ous in the S100 simulation. It is likely that this increase 
tracks the Poisson noise in the large scale clustering of the 
sources and void distribution. For the S100 simulation, the 
spectrum shows a local maximum at scales ~ 30 h^ 1 Mpc. 
Whether this maximum is real or only the effect of the 
box size is an interesting question. Only simulations with 
larger boxes will give the answer. 
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power spectrum of the differential brightness temperature 
with 4 different assumptions used to compute the spin 
temperature (see legend and main text). The spectra are 
computed for the S100 simulation. The top panel is for 
the redshift when (x a ) = 1. and the bottom panel is for 
the redshift when (xh) =0.5. 



4.2.5. Differential brightness temperature 

In this section, we will compute directly the dimensional 
power spectrum of 5Tb. Various authors, studying the 
emission regime, compute the power spectrum of 5Tb di- 
vided by the average, redshift dependent, emission tem- 

1/2 

perature To ~ 28 (^jp) ■ When we include the absorp- 
tion regime, To has no special significance. Using the ac- 
tual average 5Tb of the box would be no better : hav- 
ing both emission and absorption, it holds no information 
about the amplitude of the fluctuations and may even be 
zero at some redshift. Using the actual rms average at each 
redshift would indeed be more relevant as it does reflect 
the amplitude of the fluctuations. However it would erase 
the information about the absolute amplitude at a given 
scale, at different redshifts. Consequently, we decided to 



plot the direct dimensional power spectrum of the differ- 
ential brightness temperature. 

Fig. [15] shows the three-dimensional isotropic power 
spectrum of 5Tb as defined by Eq. 8 for the S100 simu- 
lation, at the WFCR (top panel) and the HRR (bottom 
panel), in four different cases for the computation of the 
spin temperature. In case 1, we assume that T$ ^> Tqmb, 
in case 2 T s = Tk, in case 3 we use Eq. 9 with the box av- 
eraged value of x a , and finally in case 4, we use Eq. 9 with 
the local values of x a . As expected, the larger differences 
appear at the WFCR, where the amplitude of the power 
spectrum changes by more than a factor of 10 depending 
on the case. Case 1, ignoring the possibility of absorption, 
strongly underestimates the amplitude of the power spec- 
trum, while case 2, assuming full coupling between spin 
temperature and gas kinetic temperature, overestimates 
the amplitude. The amplitudes for cases 3 and 4 are simi- 
lar (and more realistic). However, case 4 shows more power 
at small scales and less power at large scales than case 
3. Indeed, as illustrated in Fig. [131 case 3 uses a larger 
value of x a far from the sources, thus produces a stronger 
absorption in large scale voids and more power at large 
scale in the spectrum than case 4. The behaviour at small 
scales arises from the presence of the Ly-cv halos around 
the sources in case 4. Semelin et al. ()2007|) have shown that 
the profile of these halos can be as steep as r~ 3 if the gas 
has an isothermal profile. Even at distances of a few co- 
moving Mpc, where the gas profile flattens, the Ly-cv halo 
profile behaves as ~ r - 2 - 3 (Semelin et al. 120071 Chuzhoy & 
Zheng |2007p . The contribution from these sharp x a fluctu- 
ations explains the increase in the case 4 power spectrum 
at small scales. At the HRR (bottom panel) the situation 
is simpler. The Wouthuysen-Field coupling is very strong 
and case 2, 3 and 4 are almost indistinguishable. All three 
have much larger amplitude than case 1 which is related 
to the fact that absorption can produce a larger absolute 
signal than emission, thus larger fluctuations. It is interest- 
ing to notice that the difference in amplitude is especially 
large at small scales. The contribution of the regions strad- 
dling the sharp ionization front is presumably responsible 
for this feature. In case 1, the brightness temperature rises 
from ~ mK in the fully ionized region inside the ioniza- 
tion front to a moderate ~ 20 mK emission outside of the 
front. In all other cases, the temperature first rises to a 
moderate emission in the regions just outside the ioniza- 
tion front, preheated above the CMB temperature by the 
hard part of the source spectrum, then drops to strong 
absorption farther out in the cool voids. The thickness 
of the preheated regions depends strongly on the source 
modeling : in our case it is less than 1 h^ 1 Mpc . With a 
harder spectrum these regions would be thicker, and with 
X-ray sources they would be erased and the corresponding 
feature absent in the power spectrum. 

Fig. [TBI shows the 5Tb power spectrum of both the S20 
and S100 simulations at four different stages of reioniza- 
tion (the WFCR, (x H ) = 0.1, (x H ) = 0.5, and (x H ) = 
0.9), for case 1 (the four panels on the left) and case 4 (the 
four panels on the right). Let us first comment on how the 
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Fig. 16. 3D isotropic dimensional power spectrum of the differential brightness temperature for the two box sizes. The 
stage of reionization is indicated in each panel. The four left panels assume T$ 3> Tqmb while the four right panels 
compute Ts exactly using the local values of x a . 



spectra for the two simulations connect. In most cases the 
large scales in the S20 simulation have less power than 
the same scales in the S100 simulation. This is the reason 
why EoR simulations have to use large boxes : to sam- 
ple well enough the large scale clustering properties of the 
source distribution. There is however one exception, the 
( x h) = 0.5 (HRR) in case 4, where the S20 simulation has 
more power even on its large scales. This is related to the 
fact, already discussed in section 4.2.3, that the absorp- 
tion in the voids is stronger in the S20 simulation and thus 
the fluctuations are amplified. A second point is that at 
the WFCR the case 4 spectra have a more complex signa- 
ture and a larger amplitude than the case 1 spectra which 
merely track the density fluctuations. Taking absorption 
and Ly-a coupling correctly into account is essential at 
this early redshift. A third point is that we reproduce the 
behaviour noted by Mellema et al. (|2006b|) and Lidz et al. 
(2007b), that at scales of a few comoving Mpc the am- 
plitude of the power spectrum reaches a maximum at the 
HRR, then decreases. This is true for case 1, which was 
used by Lidz et al. In case 4, with the full Ly-a modeling, 
it looks like the maximum shifts to smaller values of the 
average ionization fraction. 

5. Conclusions 

In this work, we have presented simulations of the 21 cm 
emission from the EoR. The first hydrodynamic simula- 
tions were run in cosmological boxes of sizes 20 /i -1 Mpc 
and 100 hr 1 Mpc with 2 x 256 3 particles. Then, using our 
adaptive Monte-Carlo code LICORICE, we ran, in post- 
treatment, radiative transfer simulations for the ionizing 
continuum. Finally cosmological radiative transfer simula- 
tions in the Lyman-a line were performed. The latter are 
necessary to compute the neutral hydrogen spin temper- 
ature without making any assumption about the strength 
of the coupling. This approach has not been followed be- 
fore. We model the sources of both ionizing continuum and 



Lyman-a photons as stars, using a Salpeter IMF cut off at 
100 M . The 100/i _1 Mpc simulation star formation his- 
tory is calibrated on the 20/i _1 Mpc. The photon escape 
fraction is calibrated independently in both simulations 
to produce a Thomson scattering optical depth of CMB 
photons within la of the WMAP3 value, and a complete 
reionization at redshift z ~ 6. We have not included X-ray 
sources. If we had they would heat up the IGM and reduce 
the contribution of absorption vs emission in the signal. 

We first computed the evolution of a number of box- 
averaged quantities as functions of the average ionization 
fraction. We find that, with our source modeling, the av- 
erage Lyman-a coupling coefficient x a reaches 1 for small 
values of the ionization fraction (a few 10~ 3 ). However full 
coupling (x a 3> 1) is achieved only for ionization fractions 
greater than 10 _1 . This is in line with expectations. Next 
we plotted the average spin temperature and differential 
brightness temperature. It was shown that the results are 
very different depending on how the average is computed : 
directly on the spin\brightness temperature of the parti- 
cles, or from the volume-averaged values of the gas kinetic 
temperature and coupling coefficients. The latter is often 
used, but the former is more relevant. The direct average 
tends to erase the signal as regions in absorption have a 
much stronger signal than regions in emission. Maps of 
the 21 cm emission confirm that the different assumptions 
made in computing the spin temperature lead to very dif- 
ferent levels of the signal. The Ts 3> Tqmb assumption 
is not relevant at all for our choice of source modeling. 
It would be more relevant at not-too-early redshifts if we 
included X-ray sources. As expected the full coupling as- 
sumption, x a = oo, is reasonable except in the early phase 
of reionization. In this early phase, we illustrate how com- 
puting the local values of x a produces large fluctuations in 
the brightness temperature compared to using a redshift- 
dependent uniform value. 
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Finally the 3D dimensional power spectrum of the dif- 
ferential brightness temperature was plotted, with differ- 
ent assumptions in computing the spin temperature, at 
various stages of reionization and for the two box sizes. 
The first important (and expected) result is that including 
the possibility of a signal seen in absorption adds a large 
amount of amplitude to the power spectrum. This obvi- 
ously depends on the choice of the source model. Then, in 
the early phase of reionization ((xj?) < 0.1), using the lo- 
cal value of the Lyman-a: coupling changes the spectrum 
by transferring power from large scales to small scales. 
We reproduce the behaviour shown by Mellema et al. 
(|2006b|) and Lidz et al. (|2007b|) : a maximum of the am- 
plitude at (x H ) — 0.5 for the scales around 10/i _1 Mpc. 
However, when absorption and the exact computation 
of the Lyman-a coupling are included, this maximum is 
shifted to smaller ionization fractions. 

As mentioned several times, the first modification to 
our model that will be investigated is a modification of 
the source model. Including even a limited number of X- 
ray sources such as quasars would more efficiently preheat 
the regions outside of the ionization front, reducing the 
strength of the absorption regime, or even turning it to 
emission. This will be investigated in a forthcoming paper. 

Another interesting development of this work is to in- 
clude the anisotropic effect of the peculiar velocities on 
the 21cm emission (Barkana & Loeb l2005ap . Using the 
anisotropy, it is possible to separate the contribution of 
the velocities to the 21 cm fluctuations from the other con- 
tributions and to derive constraints on the cosmological 
parameters. Since this effect is important mostly in the 
early reionization phase and on small scales, it is a logical 
follow-up of including the full modeling of the Lyman-a 
coupling. 

The implication of this work for the design of the future 
radio interferometers is simply to emphasize the potential 
for the detection of the signal during the early phase of 
reionization when the fluctuations in the signal are strong. 
Unfortunately, translating " early phase" into a numerical 
value for the redshift depends crucially on the star forma- 
tion history at these high redshifts, of which very little is 
known. From the value of the Thomson scattering optical 
depth, it is likely, however, that this early phase occurs at 
z > 12. 
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